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Understanding the mechanisms underlying ion channel function from the atomic-scale requires accurate 
ab initio modelling as well as careful experiments. Here, we present a density functional theory (DFT) 
study of the ion channel gramicidin A, whose inner pore conducts only monovalent cations and whose 
conductance has been shown to depend on the side chains of the amino acids in the channel. We 
investigate the ground-state geometry and electronic properties of the channel in vacuum, focusing on 
their dependence on the side chains of the amino acids. We find that the side chains affect the ground 
state geometry, while the electrostatic potential of the pore is independent of the side chains. This study is 
also in preparation for a full, linear scaling DFT study of gramicidin A in a lipid bilayer with surrounding 
water. We demonstrate that linear scaling DFT methods can accurately model the system with reasonable 
computational cost. Linear scaling DFT allows ab initio calculations with 10,000 to 100,000 atoms and 
beyond, and will be an important new tool for biomolecular simulations. 
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1. Introduction 

Ion channels reside in cell membranes and control the passage of ions, a vital function in regulating 
cell behaviour. The gramicidin A (gA) ion channel is one of the smallest and simplest ion channel 
systems [IHSj, and was the target of the first antibiotic used in clinical practice. It is formed from 
a stack of two short beta-helix polypeptides, as shown in Fig. IJ Each polypeptide contains fifteen 
amino acids in a sequence of alternating left-handed and right-handed types. The latter helicity 
is rarely found in natural peptides, and owing to the alternating helicity, all peptide residues face 
outward into the lipid region, forming the unique helix structure with the beta type backbone 
structure. As is seen in Fig. [T] the channel structure forms a narrow pore, through which water 
molecules or ions can permeate in single file. 

In spite of its simple structure, gA shows high selectivity towards ion permeation. It is selective 
for monovalent cations, such as H + , K + , Na + , or Cs + , showing no measurable permeability 
for anions or polyvalent cations [4j. Among many ion channels, gA was one of the first whose 
atomic resolution structure was available. Many experiments have been performed to investigate 
its transport properties and to clarify the mechanism of its high selectivity. Considering the 
existence of significant experimental information as well as its simple structure, the gA system 
is of great importance as a model membrane protein j3]. 
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In the computational studies of the gramicidin- A, two high resolution structures have been 
commonly used so far. They are labelled as 1MAG and UNO in the PDB database, respectively. 
The former is obtained from solid-state NMR for the channel in DMPC membrane|5], while 
the latter one is from solution NMR, in detergent micelles [6]. While these two structures are 
topologically similar, they exhibit some variation in helix properties and residue Tryptophan 9 
(Trp9) orientation (Fig. [2^a)). The 1MAG structure helix has R=6.5 residues per helix turn and 
d=4.95 A helix turn height (pitch) in constrast to R=6.3 and d=4.79 A of the UNO structure. 
The structural disparities could have originated from any aspect of experimental procedure, with 
little to indicate the optimal channel geometry: a classical potential MD study comparing the 
two structures in a lipid membrane environment predicted that the two backbones after dynamic 
relaxation starting from two structures become almost the same and that the Trp 9 spends 80 
% of the time in the UNO orientation and 20 % in the 1MAG orient ation[7j. Although we do 
not treat the ion permeation directly in this study, it is important to note the experimental 
result that the rate of the ion permeation is significantly changed by the replacement of the 
Tryptophan(Trp) residues, which have dipole moments (~2D), with the other residues (6[ 0- 
II Oj . When it is replaced with nonpolar Phenylalanine (Phe), the ion permeation rate is greatly 
reduced. The change of the ion permeation rate is larger if the number of substituted Trp's is 
larger, and also strongly depends on which Trp's are replaced; the mutation of Trp 9 is reported 
to show the largest change. 

The system is an ideal target for theoretical studies. Although most theoretical studies were 
based on semi-microscopic models in the early stage, it is necessary to employ atomic-scale 
simulations, such as molecular dynamics (MD), for quantitative and detailed understanding of the 
permeation properties [llj. Many MD simulations of the gA ion channel system have been already 
reported f71 ll2f[T6] and nowadays it is common to treat the gA system in the channel environment, 
that is gA with surrounding lipid bilayers (a simple model for the membrane) and bulk water. 
Such MD studies have clarified many aspects of the system, but they have a serious problem 
that the results sometimes depend on the type of the force field used in the calculations [15J. 
Since the gA system with its environment requires accurate modelling of many different types 
of interaction, the transferability of classical force fields may lead to poor description of specific 
interactions. For example, water molecules and ions are transported through the narrow pore 
of the gA channel. It is unreliable to use an unpolarisable force field for water molecules or 
ions in such a unique constrained space [TT]. One of the well-known problem of the theoretical 
studies using the classical force fields is that the calculated energy barriers reported for the 
ion translocation in the gA pore are rather too high to explain the experimentally observed 
permeation rate|15 1 [TE [ [T9j. It is therefore important to apply more accurate methods to the gA 
system. 

Ab initio approaches based on density functional theory (DFT) have been playing important 
roles in clarifying the physical properties of various kinds of materials, though standard 
implementations of DFT are limited in the system size which they can model. A simplified 
gA model system has been modelled by DFT calculations to study the proton diffusion in the 
pore of the gA channel [20], A recent quantum mecahnics / molecular mechanics (QM/MM) 
study on the gA system showed the importance of polarisation of the water molecules in the 
gA pore|17| |2"T] , However, theoretical studies based on DFT are still very limited because of the 
demanding computational time and scaling with system size. Recent work has developed DFT 
methods which scale linearly with system size[22j, allowing DFT calculations on systems with 
10,000-100,000 atoms or more[23j. This brings the possibility to perform DFT calculations on 
entire biological molecules and their environment. In this work, we calculate the atomic and 
electronic structures of the isolated gA using the linear scaling DFT code CONQUEST|24[ 125) . 
In particular, we investigate the effects of the side chains of the channel on its structural and 
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Figure 1: (a)Structure of gramicidin-A ion channel in lipid bilayers. The upper and lower area are bulk water 
region, (b) Top (upper) and side (lower) views of the gramicidin-A molecule. (Online version in colour.) 



electronic properties. Although we only treat the isolated gA here, we still obtain important 
information about the structural stability and electronic structure of the gA molecule itself. We 
expect that such information would also help to clarify and improve the accuracy of a force field 
in the future. To our knowledge, it is the first DFT study to treat the whole gA molecule without 
any simplification. We will report linear scaling DFT calculations on the gA channel with full 
environment in a future publication. 

The rest of the paper is organized as follows. In Sec. [2] we briefly explain the computational 
method used in this work. The results of our calculations on the atomic and electronic structures 
of the isolated gA molecule are presented in Sec. [3] Finally, concluding remarks are given in 
Sec.H 



2. Computational Method: Local orbital and order-N methods 

The calculations in this work are performed using the CONQUEST code, which is designed for very 
large-scale DFT calculations [23 -26J. The method used in the code is based on density functional 
theory (DFT), with the pseudopotential technique. In this work, we employ the generalized 
gradient approximation (GGA) by Perdew, Burke and Ernzerhof (PBE)[27] for the exchange- 
correlation energy. 

In CONQUEST we work with local orbitals { <pi a (r) }, which are centered on the position of 
each atom i with the index of orbital a. We represent the Kohn-Sham orbitals ^ u (r) [y: band 
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index) as linear-combinations of the local orbitals. 

^^r) = ^c ia {v)(i)i <x {r) (2.1) 

ia 

These local orbitals, which we call support functions, are themselves expressed by basis set 
functions. CONQUEST provides two types of basis set functions, B-splines on regular grids for 
accurate calculations since the basis set is systematically improvable [28J, and pseudo atomic 
orbital (PAO) basis set for efficient calculations [29]. Users can choose one of the two basis set 
functions depending on the purpose. We only use the latter in this study. 

Although Conquest can solve for the Kohn-Sham orbitals using a conventional 0(N 3 ) 
diagonalisation technique, the big advantage of using CONQUEST is that it can also employ 
O(N) calculations. In the 0(N) calculations, the code uses the density matrix minimization 
method proposed by Li, Nunes and Vanderbilt (LNV)[30, 31"] . Since the detailed explanation of 
the method has already been given in our previous papers, only the main points are explained 
in the following. In the density matrix minimization method, we calculate the density matrix 
instead of the Kohn-Sham orbitals. Formally, the density matrix is defined as 

p(r,r') = Y,fuK(r)^u(r'). (2.2) 

The total energy within DFT can be expressed by the density matrix p and we optimise it 
instead of solving Kohn-Sham orbitals. In the LNV method, the matrix K is expressed as K = 
3LSL — 2LSLSL to keep the density matrix weakly idempotent. Then, we introduce a cufoff 
radius Rl for the off-diagonal elements of the L matrix in order to utilize the locality of density 
matrix. One important advantage of the LNV method is that it is variational with respect to 
Rl. We can examine the cutoff energy dependence of total energy easily and can evaluate an 
error from using a finite cutoff of Rl. This convergence behavior is reported in Sec. [3]jc]) for the 
DFT calculation of the isolated gA system. 

In both diagonalisation and 0(N) methods, the code can perform non-selfconsistent (NSC) 
DFT calculations using the Harris-Foulkes energy functional with a charge density constructed 
from the superposition of atomic charge|32| I33j. Although the NSC method is less accurate than 
usual DFT calculations with the self-consistent field (SCF) method, it enables us to perform 
efficient structure optimization or molecular dynamics. Note that the code can calculate forces 
which are completely consistent with the NSC method |2"6] I34j . This method may be useful for 
our future study on gA systems, and its accuracy is also checked in this paper. 



3. Results and Discussion 

(a) Structural Properties 

In this section, we report our results on the structural properties of the gA molecule in the 
vacuum. Here, we perform structure optimisation using two different high-resolution structures 
as initial atomic positions, 1MAG and UNO structures (Fig. [2^a)). In order to employ efficient 
structure optimisation, we first use the NSC method, then switch to the usual SCF method[34j. 
In both NSC and SCF methods, all atomic positions are relaxed using the standard conjugate 
gradient method until the maximum force component becomes smaller than 0.05 eV/ A. Note 
that the calculations in Sec. |^] and Sec. 3^ are done by the standard diagonalisation method 



with the DZP basis set, explained in Sec. |2j and GGA-PBE functional |27j. We use only a T point 



Figure 2: (a) Comparison of 1MAG (blue/dark) and UNO (red/light) experimental structures. Ab initio 
structural optimisation initial (black/dark) and final (red/light) structures in the (b) fMAG case; 
(c) f JNO case. (Online version in colour.) 



for /c-point sampling and the cutoff energy for the charge density grid is 150 Hartree. Periodic 
boundary is used with a unit cell of 23.3Ax23.3Ax35.0A. 

Figure [2^b) and (c) show the initial and optimised atomic positions starting from 1MAG 
and UNO structures, respectively. Although the original structures are obtained in different 
environment (in DMPC bilayers or in SDS micelles), we found that the optimised positions are 
close to the original ones in both cases. These figures also show that 1MAG has larger differences 
between the initial and final coordinates, especially for the backbone structure. In order to 
investigate the backbone structure in detail, we calculate the number of residues per helix turn R 
and the helix pitch d following the method in Ref. [35]. The root mean square deviations (RMSD) 
for the atoms of the backbone are also calculated. These structural parameters of 1MAG and 
UNO, before and after the structure optimisation, are listed in Table [T] The results show that 
the optimised 1MAG and UNO have different backbone structures. In addition, the change of 
the backbone structure is small in UNO case, while it is significant in 1MAG case, especially for 
helix pitch d. Although the original value of d in 1MAG case is already larger than that of UNO, 
the increase of d in 1MAG case is 0.15 A, which is larger than 0.08 A in UNO case. We observe 
a similar behaviour in the hydrogen bond lengths. The f3 -helix structure is stabilised by the 
hydrogen bonds formed between the backbone amide hydrogen and the carbonyl oxygen atoms 
parallel to the pore axis. The increase of the hydrogen bonds by the structure optimisation is 
12% in 1MAG, while it is 3% in UNO case. 

For the 1MAG structure, it is useful to know what causes the large change of the backbone 
structure during the optimisation: either the effect of backbone itself or that of side chains. 
One of the main differences between 1MAG and UNO structures is the different orientation of 
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Table 1: RMSD between the initial (I) and optimised final (F) channnel geometries and corresponding helix 
properties (number of residues per helix turn R and helix pitch d). Self consistent (upper) and non 
selfconsistent (lower) results are shown. 

the indole ring of one Tryptophan (Trp9). The atomic positions of other side chains are also 
slightly different. These differences may cause the present different backbone structures. In order 
to investigate the effects of side chains on the backbone structure, we prepare a fictitious channel 
by replacing all amino acides of the gA with alanine, keeping the same backbone structure as 
the original 1MAG structure, shown in Fig. [3] Hereafter, we refer to this system as an alanine 
tube. We have relaxed the atomic positions of this system and analyse the values R and d of 
the resulting optimised structure. As is shown in Table [TJ the change of the structure by the 
relaxation is found to be very small. The RMSD value is 0.21 A and the change of d is only 
0.04 A. This result suggests that the large changes in the original 1MAG structure is due to 
the side chains, most likely steric hindrance effects. This explains the increase in the helix pitch 
during structure optimisation. 

From the results so far, we can conclude that the original UNO structure is close to the 
optimised geometry in the gas phase, while 1MAG needs a large relaxation to reach the stable 
structure. By comparing the total energy, we find that initial UNO structure is more stable 
than 1MAG by 10.3 eV. The corresponding energy difference calculated by classical force fields is 
10.4 eV by CHARMM (v27) and 19.7 eV by AMBER (v99)[36]. The value found using AMBER 
is larger than the DFT value, but these values are essentially consistent with the DFT result. 
Following structure optimisation with DFT, this large energy difference between the 1MAG and 
UNO structures is greatly reduced, but the final structure of UNO is still more stable than that 
of 1MAG by 1.8 eV. Judging from this result in the gas phase, UNO structure seems to be more 
appropriate for the stable geometry in a lipid bilayer. Since 1MAG geometry was measured in 
ordered bilayers, some fluctuations in side-chain orientation may have been overemphasized in 
the measurement. However, we should also comment that the energy difference between these 
structures may well become smaller in a lipid bilayer environment because the thickness of the 
lipid bilayers is usually larger than the length of gA. The strength of the interactions between 
gA and lipid molecules should be different depending on the length of the gA molecule and the 
1MAG structure may be more stabilised by this effect than UNO, reducing the energy difference. 
We will address this in future work. 

Next we report the stabilisation energy of the two helices of gA. When the gA ion channel 
is open for ions, the upper and lower helices must be associated to form a dimer. They are also 
sometimes dissociated by the lateral movement in the membrane (lipid bilayers) and in such cases 
the channel is closed for the ion permeation. Thus the stabilisation energy of the two helices can 
be regarded as a measure for the stability of the open channel structure. In this respect, the 
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Figure 3: (a) Atomic structure of the gramicidin A (PDB code : 1MAG). (b) Alanine tube, which is made by 
replacing all the residues of gramicidin A by that of alanine. (Online version in colour.) 



stabilisation energy is important and we calculate the energy defined as 

Es — -E^dimcr 2E mono (3-1) 

Here, E mono is the total energy of the upper or lower peptides in isolation, whose atomic positions 
were relaxed. -Earner is the total energy of the associated gA, calculated with the basis set 
superposition error (BSSE) correction using the counterpoise method[37J. With the optimised 
UNO structure, E s is calculated to be 0.88 eV. This value is reasonable given the number of 
hydrogen bonds between the two helices, which have three N-H — O and three O-H — N bonds. 
Note that the exchange-correlation functional of PBE has been reported to be accurate for the 
expression of hydrogen bonds[38, [39] . In the membrane, this associated gA dimer structure will 
also be affected by the movement of surrounding lipid molecules. It will be important to compare 
this stabilisation energy with lipid-gA and lipid-lipid interactions in the future. 

Finally, we note that, while we have mainly discussed results obtained with the SCF method, 
the differences between the NSC and SCF methods are small for the optimised structures (see 
Table [T]). Since the NSC method is stable, and more efficient, it may play important roles in the 
future study on larger and more complicated full gA systems. 

(b) Electronic Structure 

In this section, we are particularly interested in the effects of the side chains of the ion channel 
on the electronic structure in the pore region of the gA. As mentioned in Sec. [T] it is known that 
the substitution of Trp significantly affects the rate of ion permeation[6j IHrllO]. It is important 
to understand whether such differences come from the different electronic structure of the side 
chains or some other factors, for example, from the structural effects of the interactions of the 
side chains with lipid molecules. We note that there is a report showing that the conductance of 
a monovalent cation is correlated with the electrostatic interaction energy between the ion and 
the dipole moment of the indole ring of the tryptophans |9]. 
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Figure 4: Electronic density of states for the alanine tube (red dotted line) and the 1MAG structure (blue solid 
line). (Online version in colour.) 



First, we show in Fig. [4] the density of states of gA of 1MAG optimised structure (blue), as 
well as that of alanine tube (red) introduced in the last section. We also calculated the DOS of 
the UNO structure, but the difference to the 1MAG DOS is small, particularly around the band 
gap, so this is not shown for clarity. We can see that the density of states of the 1MAG shows a 
HOMO-LUMO energy gap of 3.3 eV. The energy gap of alanine tube is much larger than that 
of gA, at around 4.2 eV, implying that HOMO and/or LUMO of gA mainly come from its side 
chains. The alanine tube also shows a very low DOS around 4eV below the Fermi level, which is 
not present in the gA DOS. 

We studied the effects of side chains on the pore of the gA by investigating the charge density 
and the electrostatic potential. Ion transport depends critically on the electrostatic potential 
inside the channel pore, which in turn depends strongly on the charge density. The potential also 
affects the dynamics of the water molecules arranged in the pore. If the charge density is plotted 
as a function of distance from the channel axis, i.e. the centre of the pore, the two systems are 
found to be almost identical until the distance is greater than 4.4 A, the approximate radius of 
the gA. We can thus conclude that the effects of the side chains on the charge distribution in 
the pore region are very small. 

We next investigate the electrostatic potential, calculated as a sum of the local pseudopotential 
and the Hartree potential using the self-consistent charge density. Figure ^a) shows a cross 
section of the potential perpendicular to the channel axis for the 1MAG structure, 8.8 A along 
the axis away from the centre of the gA. In the figure, areas of the lower potential are located 
near the peptide bonds which form hydrogen bonds. In contrast, the region near the centre of 
the pore shows a high but flat potential, making it accessible to positively charged particles. 
Figure pjb) shows the equivalent cross section for the alanine tube having a 1MAG backbone 
structure. There is a clear difference between Figs. |5ja) and (b) in the region of side chains, for 
example in the upper area. However, the two potentials are almost the same in the pore region. 
Figure [5jc) shows the difference of these two potential surfaces, implying that the two pores 
exhibit identical electrostatic potential within 0.01 eV. The same results are found at all points 
along the axis of the pore. We conclude that the effects of dipole moments of Trp's are well 
screened and the electronic structure of the side chains of the gA does not affect the electrostatic 
potential in its pore. 
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(c) Order-N calculations on the isolated gramicidin A 

It will be important in the future to perform all-atom DFT simulations of the gA channel in 
a lipid membrane with bulk water. To perform such large-scale DFT calculations, we will need 
to use an O(N) method. As a preparation for these future studies, we have made preliminary 
O(N) calculations on the isolated gA, in order to clarify the accuracy of the method on the gA 
systems. In particular, we can test the effect of the truncation of the density matrix since we can 
perform exact calculations (we have similarly tested convergence on a DNA ten-mer in water |39| 
which found essentially complete convergence for a range of lOA). Here, we use the NSC method 
because the convergence behaviour will be qualitatively the same for the SCF methods. For the 
basis sets, contracted DZP basis sets|29| are used. 

Figure [6] shows the total energy calculated for different values of the cutoff, Rl. We find that 
all calculations are stable in the optimisation of L matrix. Since the present system contains only 
552 atoms, we can employ the exact diagonalisation method and the total energy calculated by 
the method is also plotted by a dotted line in the figure. The inner panel of Fig. [6] shows the 
difference between the O(N) and exact diagonalisation methods, on a log scale. This shows the 
total energy by the O(N) method converges very fast towards the exact value. When the cutoff 
larger than 7.5 A is used, the energy discrepancy is less than 1 meV per atom. This cutoff is small 
enough to realize actual calculations. In our previous O(N) DFT study of a different system, Ge 
nano structures on Si(OOl) surface, we actually used 10.8 A for the structure optimisation of the 
systems containing more than 20,000 atoms [30]. 

We expect that similar accuracy in the total energy would be obtained in the calculations of 
full gA systems, because the water and lipid molecules have large HOMO-LUMO energy gap, 
resulting in fast convergence of the total energy with respect to R^. In fact, in our investigations 
of dry and hydrated DNA systems, we found that the difference caused by bulk water is extremely 
small[39j. Based on these results, we expect that the O(N) calculations of the full gA systems 
will be both highly accurate and useful for understanding structure and function. These studies 
are already under way and will be presented in a future publication. 
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Figure 6: Order-N band energy behaviour with increasing range of auxilliary density matrix L (inset: convergence 
to diagonalisation result). (Online version in colour.) 



4. Conclusions 

In this work, we have performed a DFT study of gramicidin A (gA) to characterise the channel 
in isolation and examined its electronic structure. In particular, we have investigated the effects 
of the side chains on the structure and electronic properties of the channel. 

For the structural stability of the gA in the gas phase, we have performed geometry 
optimisation starting from two different experimental atomic coordinates, 1MAG and UNO. 
We have found that the two calculations did not result in a common structure, but reinforced 
the difference in the helix properties, R and d. Although the original 1MAG structure has larger 
helix pitch than UNO, the optimisation of 1MAG structure exhibited a helix expansion. On the 
other hand, optimised UNO structure in the gas phase was found to be very close to the original 
geometry. To investigate the effects of side chains on the helix structure, we have introduced 
an artificial alanine tube system, which has a common backbone structure with 1MAG but 
all residues are replaced by that of alanine. We have found that the optimised helix pitch of 
the alanine tube is close to that of the original 1MAG structure. This result suggests that the 
original 1MAG structure in the gas phase includes steric hindrance between the side chains. 
Comparing the total energy of the optimised structures, UNO structure is more stable than 
1MAG structure by 1.8 eV. However, this energy difference may become smaller in the membrane 
environment. Since the thickness of lipid bilayers is usually larger than the length of gA and 
1MAG structure has smaller length mismatch, 1MAG structure may be more stabilized than 
UNO by the interactions between gA and lipid molecules. This aspect is interesting and should 
be clarified in the future study. In addition to the structure optimization of gA, we have also 
calculated the stabilisation energy for the upper and lower helices to make a dimer. This energy is 
important because the gA system shows a gating behaviour for the ion permeation by association 
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or dissociation of the two helices. The energy is calculated to be 0.8 eV, implying that reasonably 
strong hydrogen bonds are formed between the two helices. 

We have also investigated the electronic structures of gA and the alanine tube. The charge 
density and electrostatic potential in the pore of the two systems have been compared to clarify 
the effects of side chains. The potential is relevant for the ion transport and the dynamics of the 
water molecules arranged in a single file in the channel. We found that the potential in the pore 
is almost identical between the two systems. The results imply that the sensitivity of the channel 
function to the type of residues cannot be explained by the electrostatic interactions made by 
the dipole of Trp's. This does not support some of the proposals in previous reports. 

In the next step, it is desirable to simulate the system in the channel environment, that is gA in 
membrane with bulk water. We have already started to work on this large and complex gA system 
containing more than 16,000 atoms using an O(N) DFT method. The code CONQUEST used in 
this work has an ability to treat very large systems containing many thousands of atoms, with the 
density matrix minimisation method. As a preparation for such large-scale DFT studies, we have 
applied the method to the isolated gA in this paper. By calculating the total energy dependence 
on the cutoff of the auxiliary density matrix Rl, we have found that the total energy converges 
very fast to the exact value. The errors from using the finite cutoff of Rl become negligibly small 
if we use Rl larger than 10 A. It is encouraging and we believe that the method will become one 
of the standard theoretical techniques to study the atomic and electronic structures of complex 
biological systems. 
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